## pairwise comparisons using only white respondents
##
## August 5, 2019
##
## Kevin Quinn
## University of Michigan
##

library(MCMCpack)
set.seed(578833)

mydata <- read.csv("../ScaleRaceSpring2017clean.csv")


## subset data
mydata <- mydata[mydata$R.race == "White",]
mydata <- mydata[!is.na(mydata$R.race),]


## keep responses and photo IDs but nothing else
inds <- 1:24
keep.vars <- c(paste("X", inds, "_Q58", sep=""),
               paste("X", inds, "a.id", sep=""),
               paste("X", inds, "b.id", sep=""))

mydata.sub <- mydata[, keep.vars]



## reshape into long format
mydata.sub.long <- reshape(mydata.sub, direction="long",
                           varying=list(c(1:24), c(25:48), c(49:72)),
                           v.names=c("Y", "photoID.a", "photoID.b"),
                           ids=rownames(mydata.sub))


## keep only obs with non-missing values for MM
mydata.sub.long <- na.omit(mydata.sub.long)
mydata.sub.long <- mydata.sub.long[mydata.sub.long$Y != "",]

## convert to character variables
mydata.sub.long$Y <- as.character(mydata.sub.long$Y)
mydata.sub.long$photoID.a <- as.character(mydata.sub.long$photoID.a)
mydata.sub.long$photoID.b <- as.character(mydata.sub.long$photoID.b)

## recode Y to match what MCMCpaircompare expects
for (i in 1:nrow(mydata.sub.long)){
    if (mydata.sub.long$Y[i] == "the person in photo 1"){
        mydata.sub.long$Y[i] <- mydata.sub.long$photoID.a[i]
    }
    if (mydata.sub.long$Y[i] == "the person in photo 2"){
        mydata.sub.long$Y[i] <- mydata.sub.long$photoID.b[i]
    }    
}

mydata.sub.long <- mydata.sub.long[, c(5, 3, 4, 2)]

mydata.sub.long <- mydata.sub.long[sample(1:nrow(mydata.sub.long),
                                          size=1000, replace=FALSE),]

cat("\n\nN =", nrow(mydata.sub.long), "\n\n") 


## starting values and set up
pnames <- sort(unique(c(mydata.sub.long$photoID.a, mydata.sub.long$photoID.b)))

raw.counts <- table(mydata.sub.long$Y)
zero.names <- pnames[!(pnames %in% names(raw.counts))]
if (length(zero.names > 0)){
  old.names <- names(raw.counts)
  raw.counts <- c(rep(0, length(zero.names)), raw.counts)
  names(raw.counts) <- c(zero.names, old.names)
}
raw.counts <- raw.counts[pnames]
raw.ranks <- rank(raw.counts)
theta.start <- 2*(raw.ranks / length(pnames) - 0.5)


b3.pairwise.out <- MCMCpaircompare(mydata.sub.long,
                                   theta.constraints=list(mw043.J2="-",
                                                          mb017.J2="+"),
                                   theta.start=theta.start,
                                   alpha.fixed=TRUE,
                                   burnin=50000, mcmc=1000000, thin=50,
                                   verbose=10000, seed=2202029
                         )


save(b3.pairwise.out, file="MM3.2.b3.pairwise.out.Rda")
